function st = find_eigs
addpath ./build

% x_struct = load('x_opt.mat');
% a_struct = load('a_opt.mat');

% x0 = x_struct.x_opt;

pr = [1 -0.178825138445300];

x0 = [-0.399304397122687
   0.329870352753882
   0.139253019211256
  -0.316094874566801
   0.160534051402593
  -1.613330187248071
   0.386695879907189
   1.462733580986071
   2.088313203546850
  -4.615392284425383];

x0= [-0.352041903104716
   0.371883488140480
   0.043081928813040
  -0.344595567966600
   0.414925232230526
  -1.409138550715646
   0.160855859354077
   1.312163695878589
   1.701542466656633
  -6.7234832761623761];

a0=[0.599905758858752                   0                   0                   0                   0
   0.088319784190409   7.458328891437005   0.086824594333935  -2.527010921765337   0.185552626240693
  -0.182019683942709  11.618976290206323   0.150902483808764   4.214271216226255   0.316325226650619
  -0.347269593531121  -9.702266322079611  -0.134652933292173  -0.765456083469425   0.650844686308114
   0.000000000079389  10.411022920936052  -0.000000001367527   3.292172758159841   0.074934499244602];

access_A_matrix(a0);
% access_A_matrix(a_struct.a_opt);

% xr = reset_map(x0, 0);

% 
% st.jacob = jacobian(x0(2:end));

lsiObj = LSysInfo();

% CALL FSOLVE HERE
[st.x,st.fval,st.exitflag,st.output,st.jacob] = ...
    fsolve(@(x) fsl(x, x0, lsiObj), ...
    x0(2:end), optimset('Display', 'iter'));

eigs(eye(9) + st.jacob)

tth1 = get_th1(st.x, x0);
st.x = [tth1; st.x];


    function [J]=jacobian(x)
        % computes the Jacobian of a function
        n=length(x);
        fx=fsl(x);
        eps=1.e-7; % could be made better
        xperturb=x;
        for i=1:n
            xperturb(i)=xperturb(i)+eps;
            J(:,i)=(fsl(xperturb)-fx)/eps;
            xperturb(i)=x(i);
        end
    end

    function ret = fsl(x, x0, lsiObj)
        
        % Form the initial condition (on the guard)
        th1 = get_th1(x, x0);
        
        xf   = [th1; x];
        
        xr   = reset_map(xf, 0);
        data = hybrid_flow(xr, 0, 4, lsiObj);
        ret  = data.integral_curves{end}(2:end, end) - x;
    end

    function th1_ret = get_th1(xs, x0)
        [th1_ret, fval, exitflag] = ...
            fsolve(@(th1_guess) h_sca([th1_guess; xs]), ...
            x0(1), ......
            optimset('Display', 'none'));
        if exitflag ~= 1
            disp(['Error: Solving for th1 returned exit flag ' ...
                num2str(exitflag) '\n' ...
                'Cost function value: ' ...
                num2str(fval)]);
        end
    end
end
